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ABSTRACT 

Using a local model Gammie (2001) has shown that accretion discs with cooling times 
^cooi < 312“^ fragment into gravitationally bound objects, while those with cool¬ 
ing times tcooi > 311“^ evolve into a quasi-steady state. We use three-dimensional 
smoothed particle hydrodynamic simulations of protoplanetary accretion discs to test 
if the local results hold globally. We find that for disc masses appropriate for T Tauri 
discs, the fragmentation boundary still occurs at a cooling time close to Gooi = 311“^. 
For more massive discs, which are likely to be present at an earlier stage of the star 
formation process, fragmentation occurs for longer cooling times, but still within a fac¬ 
tor of two of that predicted using a local model. These results have implications not 
only for planet formation in protoplanetary discs and star formation in AGN discs, 
but also for the redistribution of angular momentum which could be driven by the 
presence of relatively massive objects within the accretion disc. 

Key words: accretion, accretion discs — planetary systems: protoplanetary discs 
— planetary systems: formation — stars: formation — stars: pre-main sequence — 
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1 INTRODUCTION 

The formation of substellar companions within protoplan¬ 
etary discs has received a great deal of attention since the 
first observation of an extra-solar planet (Mayor & Queloz 
1995). To date many additional extrasolar planets have been 
observed (see Marcy & Bntler 2000 for a review), further fu¬ 
eling the interest in the formation of planets and the evolu¬ 
tion of protoplanetary discs. The most widely studied planet 
formation mechanism is core accretion (see Lissauer 1993) 
in which planetismals grow by direct collisions to form a 
core which, when sufficiently massive (m ~ lOrriearth) then 
accretes an envelope of gas from the disc. 

An alternative mechanism for giant planet formation 
is via the gravitational instability (Kuiper 1951; Boss 1998, 
2000). This differs from core accretion in that a rocky core is 
not initially required and the process is extremely rapid. A 
Keplerian accretion disc with sound speed Cs, surface den¬ 
sity E, and epicyclic frequency n will become gravitationally 
unstable if the Toomre (1964) Q parameter 


is of order unity. A gravitationally unstable disc can either 
fragment into one or more gravitationally bound objects, 


or it can evolve into a quasi-stable state in which gravita¬ 
tional instabilities lead to the outward transport of angu¬ 
lar momentum. The exact outcome depends on the rate at 
which the disc heats up (through the dissipation of turbu¬ 
lence and gravitational instabilities) and the rate at which 
the disc cools. It has been suggested (Goldreich & Lynden- 
Bell 1965) that a feedback loop may exist where when Q is 
large cooling dominates and the disc is cooled towards insta¬ 
bility. When Q becomes sufficiently small, heating through 
viscous dissipation dominates and the disc is returned to a 
state of marginal stability. In this way Q is maintained at a 
value of ~ 1. 

Gammie (2001) has, however, shown using a local model 
that a quasi-stable state can only be maintained if the cool¬ 
ing time tcooi > where Q is the local angular frequency. 

For shorter cooling times, the disc fragments. This finding 
is consistent with Pickett et al. (1998, 2000) that ‘almost 
isothermal’ conditions are necessary for fragmentation, and 
defines a robust lower limit to the critical cooling time below 
which fragmentation occurs. It has been suggested, however, 
that self-gravitating discs strictly require a global treatment 
(Balbus & Papaloizou 1999), and while global effects are 
highly unlikely to stabilize a locally unstable disc, they could 
well allow fragmentation within discs that would locally be 
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stable. In this paper we use a global model to test whether 
a gaseous accretion disc still fragments for tcooi < 
and attains a quasi-stable state for tcooi > Although 

our choice of parameters is more appropriate for protoplan¬ 
etary accretion discs, the results may also be applicable for 
star formation in AGN discs, which are expected to become 
gravitationally unstable at large radii (Shlosman & Begel- 
man 1989; Goodman 2002), and for star formation in the 
gaseous regions of spiral galaxies (Kim & Ostriker 2002). 


2 NUMERICAL SIMULATIONS 

2.1 Smooth particle hydrodynamics code 

The three-dimensional simulations presented here were per¬ 
formed using smoothed particle hydrodynamics (SPH), a 
Lagrangian hydrodynamics code (e.g., Benz 1990; Mon¬ 
aghan 1992). In this simulation the central star is modelled 
as a point mass onto which gas particles can accrete if they 
approach to within the sink radius, while the gaseous disc is 
simulated using 250,000 SPH particles. Both point masses 
and gas use a tree to determine neighbours and to calculate 
gravitational forces (Benz et al. 1990), and the point mass 
representing the central star is free to move under the grav¬ 
itational influence of the disc gas. Although the code can in 
principle continue running if the disc starts to fragment and 
high density regions form, it quickly becomes too slow to re¬ 
alistically continue in this manner. To continue simulating a 
fragmenting disc, point masses are created (Bate et al. 1995) 
if the high density regions are gravitationally bound (gravi¬ 
tational potential energy at least twice the thermal energy). 
As with the point mass representing the central star, point 
masses within the disc may continue to accrete gas particles 
if they fall within the sink radius. An additional saving in 
computational time is also made by using individual, par¬ 
ticle time-steps (Bate et al. 1995; Navarro & White 1993). 
The time-steps for each particle are limited by the Gourant 
condition and a force condition (Monaghan 1992). 

2.2 Initial conditions 

We consider a system comprising a central star, modelled 
as a point mass with mass M*, surrounded by a gaseous 
circumstellar disc with mass Mdiac- Most of the simulations 
presented here considered a disc mass of Maisc = O.IM*, 
but a few were performed using Maisc = 0.25M*. The disc 
temperature is taken to have an initial radial profile of T oc 

® (e.g., Yorke & Bodenheimer 1999) and the Toomre Q 
parameter is assumed to be initially constant with a value 
of 2. A stable accretion disc where self-gravity leads to the 
steady outward transportation of angular momentum should 
have a near constant Q throughout. A constant Q together 
with equation ( 1 ) then gives a surface density profile of E oc 
and hydrostatic equilibrium in the vertical direction 
gives a central density profile of p oc r~^. 

The disc is modelled using 250,000 SPH particles, which 
are initially randomly distributed in such a way as to give 
the specified density profile between inner and outer radii of 
Tin and rout respectively. 

The calculations performed here are essentially scale 
free. In code units, we take M* = 1, rin = 1 and rout = 25. 


If we were to assume a physical mass scale of IMq and a 
length scale of 1 an, the central star would have a mass of 
IM 0 , the circumstellar disc would have a mass of O.IMq, or 
O. 25 M 0 , and would extend from 1 au to 25 an, and 1 year 
would equal 27r code units. 

2.3 Cooling 

Since the aim of this work is to test whether the results 
obtained by Gammie (2001) using a local model still hold 
globally, we use the same approach in our global model as 
used in the local model. We use an adiabatic equation of 
state, with adiabatic index 7 = 5/3, and allow the disc gas 
to heat up due to both PdV work and viscous dissipation. 
Cooling is implemented by adding a simple cooling term to 
the energy equation. Specifically, for a particle with internal 
energy per unit mass Ui, 

— = -R!±- ( 2 ) 

dt tcooi 

where, as in Gammie (2001), tcooi is given by (3Qr^ with the 
value of (5 varied for each run. 

Although the use of the above cooling time is essen¬ 
tially chosen to compare the local model results with results 
using a global model, it can also be related (at least ap¬ 
proximately) to the real physics of an accretion disc. For an 
optically thick disc in equilibrium, the cooling time is the ra¬ 
tio of the thermal energy per unit area to the radiative losses 
per unit area. It can be shown (e.g., Pringle 1981) that in 
such a viscous accretion disc, the cooling time is given by 

= 97(7 - 1 ) ^ 

where 7 is the adiabatic index, II is the angular frequency, 
and a is the Shakura & Sunyaev (1973) viscosity parameter. 

3 RESULTS 

3.1 Mdisk = O.IM. 

We use a global model to consider how cooling times of 
tcooi = and tcooi = 312”^ affect the gravitational sta¬ 

bility of a viscous accretion disc with a mass of Maisc = 
O.IM*. Figure ^ shows the disc surface density, E, at the 
beginning and end of the tcooi = 512”^ simulation. Since we 
have not attempted to model the inner boundary condition 
in any detail, there is rapid accretion and a drop in surface 
density close to the inner boundary. Apart from particles 
with radii between 1 and 2 being accreted onto the central 
star, the surface density profile does not change significantly 
during the course of the simulation. This also occurs for 
tcooi = 312“^ and a corresponding figure is consequently not 
shown. 

Figure ^ shows the final equatorial density structure of 
the tcooi = simulation. The central star (not shown) is 

located in the middle of the figure and the x and y axes both 
run from -25 to 25. Figure]^ shows the Toomre Q parameter 
for the same simulation at the beggining (t = 0 ), approx¬ 
imately a third of the way into the simulation (t = 876), 
and at the end (t = 2932) of the simulation. After t = 2932 
time units, particles at the outer edge of the disk (r = 25) 
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Figure 1. Surface density for tcool = and Mdisc = O.IM* 

at the beginning and end of the simulation. Apart from particles 
with initial radii between 1 and 2 being accreted onto the central 
star, the surface density does not change significantly during the 
simulation. 



Figure 2. Equatorial density structure for tcool = and 

A^disc = O.IM*. The disc is highly structure with the instabil¬ 
ity existing at all radii. The density has, however, not increased 
significantly and the disc is in a quasi-stable state with heating 
through viscous dissipation balancing cooling. 


have completed 3.7 orbits while particles at r = 1 have com¬ 
pleted 466 orbits. The disc is highly structured and the in¬ 
stability exists at all radii. However, at no point in the disc 
has the density increased significantly and no fragmentation 
has taken place. Figure ^ shows that at t = 2932 (dashed 
line) the Toomre Q parameter, which initially had a con¬ 
stant value of 2 (solid line), is close to 1 for radii between 1 
and 15. Comparing Q at t = 876 and at t = 2932 suggests 
that the cooling may initially reduce Q to below 1, causing 
the gravitational instability to grow, heating the disc and 
returning Q to a value of order unity. Figure ^ also shows 
this region of low Q moving to larger radii with time. This 
is a consequence of both the cooling time and the dynamical 
time (the timescale on which heating occurs) both depend¬ 
ing directly on radius. The instability therefore starts at the 
inner radii and by the end of the simulation has reached 
the edge (r = 25) of the disc. The disc has reached a quasi¬ 
steady state in which cooling is balanced by heating through 
viscous dissipation resulting from the growth of the gravita¬ 
tional instability. 

Figure ^ shows the final equatorial density structure of 
the tcool = 317“'^ simulation. We were only able to run this 
simulation for a total of 504 time units and hence Figure 
^ shows only the inner 8 radii of the simulation. The cen¬ 
tral star is again in the middle of the figure. Figure shows 
the Toomre Q parameter at three different times during the 
simulation. The initial Q for tcooi = 3f7~^ is the same as 
for tcool = 517 and is therefore not shown in Figure ^ 
Figure ^ shows that not only is the disc highly structured, 
the bright dots also indicate that fragmentation is taking 
place. The fragments that can be seen in Figure ^ are all 
gravitationally bound. To reach the time shown in Figure ^ 
it was necessary to convert most of the fragments into point 
masses (Bate et al. 1995). Figure ^ shows that at t = 192 
Q is less than 1 for radii between 1 and 5. At later times 



Figure 3. Toomre Q parameter at the beginning (t = 0), one 
third of the way (t = 876), and at the end {t = 2932) of the tcool = 
512“^, Mjisc = O.IM* simulation. At the end of the simulation 
Q ~ 1 at radii between 1 and 15. A low Q region also moves to 
larger radii with time, indicating the current radius of maximum 
instability growth. 

the minimum Q value moves to larger radii and the value 
of Q at smaller radii increases to a value of between 1 and 
1.5. As in the tcooi = 517“^ simulation, the gravitational 
instability grows from the inner radii and moves to larger 
radii with time. While Q is below 1, the instability grows 
rapidly, producing gravitationally bound fragments. These 
fragments tidally interact with the disc, heating it and caus¬ 
ing Q to increase to a value greater than unity. Once Q > 1 
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Figure 4. Equatorial density structure for tcool = 3^2“^ and 
A^disc = O.IM*. The disc is highly unstable and is fragmenting. 
The fragments are all gravitationally bound. 


LO 



Figure 5. Toomre Q parameter at times of t = 102, t = 362, and 
t = 504 for fcool = 312“^ and Mmac = 0.1M». The value of Q is 
reduced such that fragmentation of the disc takes place. A region 
of low Q moves to larger radii with time, indicating the current 
fragmentation radius. 

fragmentation ceases, and the disc, at that radius, becomes 
gravitationally stable. There is, however, likely to be a sig¬ 
nificant amount of angular momentum transport driven by 
tidal interactions between fragments and the gaseous disc 
(Larson 2002). 

Simulations using tcool = 10f2“^ and tcooi = were 

also performed. Both simulations were run for more than 
2000 time units, almost 4 times longer than the tcooi = 
simulation. The tcool = 1017“^ simulation showed no notice¬ 
able structure and the Toomre Q settled down to a value 


between 1 and 2. For tcool = 417“^ the disc structure was 
similar to that obtained using tcool = 5f7“^ except for the 
presence of a few high density clumps at radii of ~ 10 and 
greater. By the end of the simulation (t = 2136) one of 
the clumps, at a radius of 9.5, had satisfied the condition 
for point mass creation, i.e. it was gravitationally bound. 
This would suggest that, globally, fragmentation may occur 
at slightly greater cooling time than expected using a local 
model (Gammie 2001). However, unlike the tcool = 317“'^ 
simulation, the tcool = 417“^ simulation did not fragment at 
all radii, with the inner ~ 10 radii of the disc remaining in 
a quasi-stable state. 

3.2 Mdisc = 0.25M. 

A simulation using a heavier disc (Mdiac = 0.25M*) was 
also performed. Although most T Tauri discs have masses 
significantly less than this (Beckwith et al. 1990), there are a 
few with comparable masses. Such a simulation may also be 
appropriate for an earlier stage of the star formation process 
when the discs are expected to be heavier. We considered 
cooling times of tcool = 1017“^ and tcool = 517“^. The tcool = 
1017“^ simulation was run for 1024 time units, at the end 
of which, as for Mdisc = O.IM*, the disc showed minimal 
structure and there was no sign of any density enhancement 
or fragmentation. 

The tcool = 517“'^ simulation was run for 684 time units 
and the final equatorial density structure is shown in Fig¬ 
ure ^ The disc is highly structured and the spirals, consis¬ 
tent with the increased disc mass (Nelson et al. 1998), are 
somewhat less filamentary then for the Mdisc = O.IM* sim¬ 
ulation with the same cooling time. Unlike the equivalent 
Mdisc = O.IM* simulation (see Figure ^) there are, how¬ 
ever, a number of high density clumps present in the disc. 
The routine for converting these clumps into point masses 
checks the nearest oo 50 SPH particles to determine if they 
are gravitationally bound (see Bate et al. 199^. Using this 
technique, the high density regions in Figure M were found 
to be gravitationally unbound and could not be converted 
into point masses. The simulation was eventually stopped 
at t = 684 when the maximum density had increased to 
a value 10^ times greater than the initial maximum den¬ 
sity. Restricting the determination of the boundness of the 
clump to only the nearest oo 50 SPH particles is somewhat 
arbitrary. By including the nearest oo 250 SPH particles, 
rather than only the nearest ~ 50, the densest clump was 
found to be just gravitationally bound, i.e. the gravitational 
potential energy was almost exactly twice the thermal en¬ 
ergy. This would imply that in more massive discs, global 
effects could act to make the disc more unstable, allowing 
gravitationally bound fragments to grow for cooling times 
greater than that obtained using a local model. This may 
have implications both for the formation of planets, or bi¬ 
nary companions (Adams, Ruden & Shu 1989), reasonably 
early in the star formation process. 

3.3 Cooling time and effective viscosity 

By using a radially dependent cooling time in our simula¬ 
tions, Equation (^) implies that the discs should settle into 
a state in which the effective viscous a (Shakura & Sunyaev 


© 0000 RAS, MNRAS 000, 000-000 




















Cooling in self-gravitating protoplanetary discs 5 



Figure 6. Equatorial density structure for tcool = and for 

a disc mass = 0.25M*. There are signs of fragmentation 

with the most massive fragment being gravitationally bound. 

1973) depends inversely on this cooling time. Since none of 
our simulations were run to a steady state, which would 
take several viscous times to achieve, it is not straightfor¬ 
ward to determine whether the local scaling of a with tcooi 
is obeyed. This is further complicated in the tcool = 

Mdisc = O.IM* simulation since the collapsed mass frac¬ 
tion increases throughout the run. We are, however able 
to compare the tcooi = 511“^, Maisc = O.IM* simulation 
with the tcool = Maisc = O.IM* simulation. Figure 

Q shows the magnitude of the mass transfer rate through 
a radius r = 15, plotted against time in code units, for 
tcool = (solid line) and tcool = 10 f2~^ (dashed line) 

and a disc mass, in both cases, of Maisc = O.IM*. Figure 
^ only considers the mass transfer rate once the instability 
has almost saturated at the radius considered. Since neither 
simulation has reached a steady state and since both discs 
are self-gravitating, the mass transfer rate can be negative 
or positive. Figure W is consequently a 60 time unit running 
average of the magnitude of the mass transfer rate and shows 
that the mass transfer rate for tcool = is, as expected, 

higher than that for tcool = 1017“^. By the end of simula¬ 
tion the mass transfer rates are both reasonably constant 
and differ by a factor of between 2 and 3. Although this 
difference is consistent with that expected from Equation 
(^), the complications in trying to accurately determine the 
relationship between a and tcool in these global simulations 
leads us to simply conclude that, as expected, the effective 
viscosity does indeed increase with decreasing cooling time. 


4 CONCLUSION 

Using a local model, Gammie (2001) has shown that for 
cooling times tcool < 3U“^ a disc will fragment into one or 
more gravitationally bound objects, while for longer cooling 
times the disc will settle into a quasi-stable state with heat¬ 
ing through viscous dissipation balancing cooling. We have 



time 


Figure 7. A running average of the magnitude of the mass trans¬ 
fer rate for tcool = and tcool = and for a disc mass 

Mdisc = O.IM*. 


performed three-dimensional disc calculations using SPH to 
test if the above still holds globally. We use the same cooling 
function as used by Gammie (2001) which, although fairly 
simplistic, can also be justified physically (Pringle 1981). 

For a disc mass of Maisc = O.IM* and for cooling times 
of tcool = 517“^ and tcool = 1017“^ we find that the disc 
settles into a quasi-stable state. At the end of the tcool = 
517“^ run, the Toomre Q parameter is of order unity between 
radii of 1 and 15 and the disc is highly structured. The 
density is, however, not significantly greater than the initial 
density and no fragmentation has occured. A comparison 
between the tcool = 5f7“^ and the tcool = 10f7“^ simulations 
also shows that once the instability has saturated at a given 
radius, the magnitude of the mass transfer rate was greatest 
for the smaller cooling time. Although we were unable to 
quantify the relationship between the cooling time and the 
viscous Q (Shakura & Sunyaev 1973), this result shows that, 
as expected, the effective viscosity does indeed increase with 
decreasing cooling time. 

For tcool = 317“^ the disc rapidly becomes unstable and 
produces numerous gravitationally bound fragments. The 
growth of these fragments starts near the inner boundary of 
the disc and had we been able to continue with the simu¬ 
lation, it seems likely that fragmentation would have con¬ 
tinued throughout the disc. A tcool = 417“^ simulation also 
produced gravitationally bound fragments but these were 
located at radii of ~ 10 or greater. The inner radii of the 
disc remained in a quasi-steady state. 

A simulation using a heavier disc (Maisc = 0.25M*) 
was found to fragment for cooling times of tcool = 517“^, 
suggesting that as the disc mass increases, global effects 
may act to make the disc more unstable. A cooling time of 
fcooi = 1017“^, however, showed no signs of fragmentation 
and hence the fragmentation boundary, even for the heavier 
disc, is still within a factor of ~ 2 of that obtained using a 
local model (Gammie 2001). For T Tauri discs, which gen¬ 
erally have masses Mdisc < O.IM* (Beckwith et al. 1990), 
the fragmentation boundary is therefore likely to be close to 


© 0000 RAS, MNRAS 000, 000-000 





6 Rice et al. 


that determined using the local model. Earlier in the star 
formation process, when the discs are expected to be more 
massive, fragmentation may occur for cooling times some¬ 
what longer than that predicted by the local model results 
(Gammie 2001). 

For quiescent T Tauri discs, a is conventionally esti¬ 
mated to be of the order of 10~^ or smaller (Hartmann et 
al. 1998; Bell & Lin 1994). The simple estimates quoted 
earlier would suggest that such discs are comfortably sta¬ 
ble. Boss (2001), however, has shown using an approximate 
treatment of the disc heating and cooling that there may 
be periods when the cooling time is comparable to the or¬ 
bital period. Our results suggest that if the disc is fairly 
massive, such short cooling times open a window of oppor¬ 
tunity for the formation of substellar objects, probably in 
the form of a multiple system (Armitage & Hansen 1999). 
Not only does this have implications for planet formation 
in protoplanetary discs, it also has implications for angular 
momentum transport. The presence of a number of massive 
clumps within a protoplanetary disc is likely to significantly 
enhance the transport of angular momentum through tidal 
interactions (Larson 2002). Similar considerations apply to 
AGN discs. Numerous theoretical studies (e.g. Glarke 1988; 
Kumar 1999; Hure 2000; Menou & Quataert 2001) show 
that AGN discs invariably become self-gravitating at radii 
of ~ 0.1 pc. Fragmentation would lead to the formation of 
stars (Gollin & Zahn 1999), and truncation of the inner ac¬ 
cretion disc (Goodman 2002). Instability of gaseous discs at 
much larger radii of 10 — 1000 pc may also occur (Shlosman 
& Begelman 1989). Fragmentation of these larger scale discs 
would yield a flattened stellar system, while rapid angular 
momentum transport could play a role in replenishing the 
inner galactic resevoir (Shlosman & Begelman 1989). 
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